パッケージとデータの読み込み

# pacmanパッケージを利用しない場合
# library(tidyverse)
# library(MatchIt)
# ...

pacman::p_load(tidyverse, 
               MatchIt, 
               WeightIt, 
               BalanceR, 
               cobalt, 
               summarytools,
               modelsummary,
               kableExtra)

# cobaltパッケージが提供するデータセットの読み込み
data("lalonde", package = "cobalt")

lalonde_df <- lalonde

マッチングにおける古典的なデータセット、lalondeを読み込みました。data(lalonde, package = "cobalt")を入力するだけで、{cobalt}パッケージ内のlaondeという名前のデータフレームが作業環境上に読み込まれます1。このデータをDay2_dfという名のオブジェクトとして改めて保存しておきます。

それでは、データの中身を確認してみましょう。普通にコンソール上でlalondeを打っても良いですが、DTパッケージのdatatable()を使うと見やすい表形式で出力されます。

# 今後、datatable()関数を使う予定はないので、DTパッケージを読み込まず、
# DT::datatable()を使用する
DT::datatable(lalonde_df)

分析に入る前に、名目変数である人種（race）をダミー変数にします。raceは3種類の値で構成されているため、生成するダミー変数も3つとなります。

lalonde_df <- lalonde_df %>%
  mutate(white    = if_else(race == "white",  1, 0),
         black    = if_else(race == "black",  1, 0),
         hispanic = if_else(race == "hispan", 1, 0),
         # 作成された以上の3変数をraceの後へ
         .after   = race)

続いて、{summarytools}記述統計量を確認してみましょう。Rコンソール上で確認するならdescr()を、R Markdown内に埋め込むならdfSummary()が良いでしょう。

# descr()を使う場合
descr(lalonde_df)
## Descriptive Statistics  
## lalonde_df  
## N: 614  
## 
##                        age    black     educ   hispanic   married   nodegree       re74       re75
## ----------------- -------- -------- -------- ---------- --------- ---------- ---------- ----------
##              Mean    27.36     0.40    10.27       0.12      0.42       0.63    4557.55    2184.94
##           Std.Dev     9.88     0.49     2.63       0.32      0.49       0.48    6477.96    3295.68
##               Min    16.00     0.00     0.00       0.00      0.00       0.00       0.00       0.00
##                Q1    20.00     0.00     9.00       0.00      0.00       0.00       0.00       0.00
##            Median    25.00     0.00    11.00       0.00      0.00       1.00    1042.33     601.55
##                Q3    32.00     1.00    12.00       0.00      1.00       1.00    7891.93    3254.81
##               Max    55.00     1.00    18.00       1.00      1.00       1.00   35040.07   25142.24
##               MAD     8.90     0.00     1.48       0.00      0.00       0.00    1545.36     891.86
##               IQR    12.00     1.00     3.00       0.00      1.00       1.00    7888.50    3248.99
##                CV     0.36     1.24     0.26       2.75      1.19       0.77       1.42       1.51
##          Skewness     1.02     0.43    -0.72       2.37      0.34      -0.54       1.58       2.29
##       SE.Skewness     0.10     0.10     0.10       0.10      0.10       0.10       0.10       0.10
##          Kurtosis     0.14    -1.82     1.62       3.64     -1.89      -1.71       1.95       7.13
##           N.Valid   614.00   614.00   614.00     614.00    614.00     614.00     614.00     614.00
##         Pct.Valid   100.00   100.00   100.00     100.00    100.00     100.00     100.00     100.00
## 
## Table: Table continues below
## 
##  
## 
##                         re78    treat    white
## ----------------- ---------- -------- --------
##              Mean    6792.83     0.30     0.49
##           Std.Dev    7470.73     0.46     0.50
##               Min       0.00     0.00     0.00
##                Q1     237.91     0.00     0.00
##            Median    4759.02     0.00     0.00
##                Q3   10923.35     1.00     1.00
##               Max   60307.93     1.00     1.00
##               MAD    7055.72     0.00     0.00
##               IQR   10655.31     1.00     1.00
##                CV       1.10     1.52     1.03
##          Skewness       1.55     0.86     0.05
##       SE.Skewness       0.10     0.10     0.10
##          Kurtosis       4.33    -1.26    -2.00
##           N.Valid     614.00   614.00   614.00
##         Pct.Valid     100.00   100.00   100.00
# dfSummary()はHTML形式で出力されるため、R Markdownに埋め込む
dfSummary(lalonde_df, style = "grid", 
          plain.ascii = FALSE, graph.magnif = 0.85) %>%
  print(method = "render", headings = FALSE)
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 treat [integer] Min : 0 Mean : 0.3 Max : 1
0:429(69.9%)
1:185(30.1%)
 614 (100.0%) 0 (0.0%)
2 age [integer] Mean (sd) : 27.4 (9.9) min < med < max: 16 < 25 < 55 IQR (CV) : 12 (0.4) 40 distinct values 614 (100.0%) 0 (0.0%)
3 educ [integer] Mean (sd) : 10.3 (2.6) min < med < max: 0 < 11 < 18 IQR (CV) : 3 (0.3) 19 distinct values 614 (100.0%) 0 (0.0%)
4 race [factor] 1. black 2. hispan 3. white
243(39.6%)
72(11.7%)
299(48.7%)
 614 (100.0%) 0 (0.0%)
5 white [numeric] Min : 0 Mean : 0.5 Max : 1
0:315(51.3%)
1:299(48.7%)
 614 (100.0%) 0 (0.0%)
6 black [numeric] Min : 0 Mean : 0.4 Max : 1
0:371(60.4%)
1:243(39.6%)
 614 (100.0%) 0 (0.0%)
7 hispanic [numeric] Min : 0 Mean : 0.1 Max : 1
0:542(88.3%)
1:72(11.7%)
 614 (100.0%) 0 (0.0%)
8 married [integer] Min : 0 Mean : 0.4 Max : 1
0:359(58.5%)
1:255(41.5%)
 614 (100.0%) 0 (0.0%)
9 nodegree [integer] Min : 0 Mean : 0.6 Max : 1
0:227(37.0%)
1:387(63.0%)
 614 (100.0%) 0 (0.0%)
10 re74 [numeric] Mean (sd) : 4557.5 (6478) min < med < max: 0 < 1042.3 < 35040.1 IQR (CV) : 7888.5 (1.4) 358 distinct values 614 (100.0%) 0 (0.0%)
11 re75 [numeric] Mean (sd) : 2184.9 (3295.7) min < med < max: 0 < 601.5 < 25142.2 IQR (CV) : 3249 (1.5) 356 distinct values 614 (100.0%) 0 (0.0%)
12 re78 [numeric] Mean (sd) : 6792.8 (7470.7) min < med < max: 0 < 4759 < 60307.9 IQR (CV) : 10655.3 (1.1) 457 distinct values 614 (100.0%) 0 (0.0%)

Generated by summarytools 0.9.9 (R version 4.1.0)
2021-07-27

まずは単純差分から

treat変数の値ごとに結果変数であるre78の平均値を計算してみましょう。

Diff_Mean_df <- lalonde_df %>% 
    group_by(treat) %>%
    summarise(Outcome = mean(re78),
              .groups = "drop")

Diff_Mean_df
## # A tibble: 2 x 2
##   treat Outcome
##   <int>   <dbl>
## 1     0   6984.
## 2     1   6349.

この結果を可視化するコードは以下の通りです。

Diff_Mean_df %>%
  ggplot() +
  geom_bar(aes(x = treat, y = Outcome), 
           stat = "identity", width = 0.5) +
  geom_label(aes(x = treat, y = Outcome,
                 label = round(Outcome, 3))) +
  labs(x = "Treatment",
       y = "Outcome (US Dollars)") +
  # scale_x_continuous()を使って0/1をControl/Treatmentに置換する
  # 目盛りはX軸上の0と1、各目盛りのラベルはControlとTreatmentに
  scale_x_continuous(breaks = c(0, 1), labels = c("Control", "Treatment")) +
  coord_cartesian(xlim = c(-0.5, 1.5))

treat == 0の回答者、つまり職業訓練を受けていない回答者の平均所得は約6984ドル、treat == 1の回答者、つまり職業訓練を受けた回答者の平均所得は約6394ドルです。その差は約-650ドルですが、職業訓練を受けた回答者の方が低所得になっています。これは直感的に納得できる結果ではないでしょう。むろん、実際、職業訓練が所得を減らす可能性もありますが、今回の結果はより詳しく分析してみる価値はあるでしょう。

ちなみに、以上の結果は単回帰分析からも確認可能です (ただし、統計的に有意ではない)。

Raw_Fit1 <- lm(re78 ~ treat, data = lalonde_df)
# modelsummaryを使った推定結果の要約
# デフォルトでは点推定値と標準誤差が出力されるが、
# ここでは点推定値と95%信頼区間を出力する。
modelsummary(Raw_Fit1,
             statistic = "[{conf.low}, {conf.high}]",
             conf_level = 0.95)
Model 1
(Intercept) 6984.170
[6275.791, 7692.549]
treat -635.026
[-1925.544, 655.492]
Num.Obs. 614
R2 0.002
R2 Adj. 0.000
AIC 12698.7
BIC 12712.0
Log.Lik. -6346.371
F 0.934

この直感的でない結果は、もしかしたらセレクションバイアスが原因かも知れません。職業訓練の対象が元々非常に所得が低い被験者になっている可能性があります。たとえば、下の図のように職業訓練の有無が教育水準や人種、これまでの所得などと関係しているとしましょう。これらの要因は回答者の現在所得にも関係していると考えられます。この場合、処置有無と所得の間には内生性が存在することになります。

本当にそうなのか、共変量のバランスチェックをしてみましょう。もし、処置有無によって回答者の社会経済的要因に大きな差があれば、内生性が存在する証拠になります。ここでは誰かが作成しました{BalanceR}パッケージを使います。

Blc_Chk <- lalonde_df %>%
  BalanceR(group = treat, cov = c(age, educ, white:re75))

{BalanceR}パッケージで共変量を指定する際、:演算子が使用可能です。white:re75は、データセットのwhiteからre75変数までをすべて指定することを意味します。Rコンソール上で変数名のみを出力するにはnames()関数を使用します。

names(lalonde_df)
##  [1] "treat"    "age"      "educ"     "race"     "white"    "black"   
##  [7] "hispanic" "married"  "nodegree" "re74"     "re75"     "re78"

whiteからre75までblack、hispanic、marriedなどの変数がありますが、これらを一々指定せず、:を使えば一気に指定することができます。それではバランスチェックの結果を確認してみましょう。

Blc_Chk
##   Covariate   Mean:0     SD:0   Mean:1     SD:1   SB:0-1
## 1       age   28.030   10.787   25.816    7.155   24.190
## 2      educ   10.235    2.855   10.346    2.011   -4.476
## 3     white    0.655    0.476    0.097    0.297  140.799
## 4     black    0.203    0.403    0.843    0.365 -167.083
## 5  hispanic    0.142    0.350    0.059    0.237   27.740
## 6   married    0.513    0.500    0.189    0.393   72.076
## 7  nodegree    0.597    0.491    0.708    0.456  -23.549
## 8      re74 5619.237 6788.751 2095.574 4886.620   59.575
## 9      re75 2466.484 3291.996 1532.055 3219.251   28.700
Covariate Mean:0 SD:0 Mean:1 SD:1 SB:0-1
age 28.030 10.787 25.816 7.155 24.190
educ 10.235 2.855 10.346 2.011 -4.476
white 0.655 0.476 0.097 0.297 140.799
black 0.203 0.403 0.843 0.365 -167.083
hispanic 0.142 0.350 0.059 0.237 27.740
married 0.513 0.500 0.189 0.393 72.076
nodegree 0.597 0.491 0.708 0.456 -23.549
re74 5619.237 6788.751 2095.574 4886.620 59.575
re75 2466.484 3291.996 1532.055 3219.251 28.700

いくつか怪しい箇所があります。たとえば、treat == 0の回答者において黒人の割合は約20%ですが、treat == 1のそれは約85%です。つまり、黒人ほどより職業訓練を受ける傾向があることを意味します。実際、標準化バイアスは-167という、非常に大きい数値を示しています。この結果を表としてまとめてみましょう。

# 標準化バイアスの可視化
# 絶対値変換。SB = 25に破線
plot(Blc_Chk, abs = TRUE, vline = 25)

かなり緩めの基準である25を採用しても、人種（黒人）、結婚有無、74年の所得のバランスが非常によくありません。内生性があると判断して良いでしょう。以下では様々な方法を使い、この内生性に対処していきたいと思います。

内生性への対処

重回帰分析

つづいて、重回帰分析もしてみましょう。用いる共変量は年齢、教育水準、黒人ダミー、ヒスパニックダミー、既婚ダミー、学位なしダミー、74・75年の所得です。lm()関数を使用し、78年の所得をこちらの変数に回帰させます。

\[ \begin{align} \textrm{re78} = & \beta_0 + \beta_1 \textrm{treat} + \beta_2 \textrm{age} + \beta_3 \textrm{educ} + \\ & \beta_4 \textrm{black} + \beta_5 \textrm{hispanic} + \beta_6 \textrm{married} + \beta_7 \textrm{nodegree} + \beta_8 \textrm{re74} + \beta_9 \textrm{re75}. \end{align} \]

Raw_Fit2 <- lm(re78 ~ treat + age + educ + black + hispanic + married + 
                   nodegree + re74 + re75, data = lalonde_df)

modelsummary(Raw_Fit2)
Model 1
(Intercept) 66.515
(2436.746)
treat 1548.244
(781.279)
age 12.978
(32.489)
educ 403.941
(158.906)
black -1240.644
(768.764)
hispanic 498.897
(941.943)
married 406.621
(695.472)
nodegree 259.817
(847.442)
re74 0.296
(0.058)
re75 0.232
(0.105)
Num.Obs. 614
R2 0.148
R2 Adj. 0.135
AIC 12617.5
BIC 12666.1
Log.Lik. -6297.752
F 11.636

共変量を統制したら処置変数の係数は約1548.244ドルです。先ほどの単回帰分析の結果とは違って、統計的に有意な正の効果が確認されていますね。

マハラノビス最近傍マッチング

重回帰分析は非常にシンプルで便利な分析方法ですが、いくつかの欠点があります。まず、重回帰分析は変数間の関係（線形結合）および誤差項の分布（平均0の正規分布）などを仮定したパラメトリック分析になります。この場合、同じ共変量を持たないケースであっても、勝手に予測を行います。重回帰分析における処置変数の解釈は「他の共変量がすべて同じ」場合の処置効果になります。これは、共変量がすべて同じ場合における（最初に見た）単純差分のようなものです。しかし、「他の共変量がすべて同じ」ケースが存在しない可能性があります。特に、共変量が多く、連続変数の場合、共変量がすべて同じことは実質あり得ないか、非常に少ないケースに限定されます。一方、マッチングを行うと、「他の共変量がすべて同じ」、または「非常に似ている」ケース間で比較を行います。

ここでは以下の4つの手法を紹介します。

  1. マッチング
    • マハラノビス距離最近傍マッチング
    • Coarsened Exact Matching (CEM)
    • 傾向スコアマッチング
  2. 逆確率重み付け

まずは、マハラノビス最近傍マッチングから始めましょう。だいたいのマッチング手法は{MatchIt}パッケージで解決できます。マッチングデータセットを作成する関数はmatchit()関数であり、使い方は

matchit(処置変数 ~ 共変量1 + ... + 共変量k, 
            data = データフレーム名, estimand = "ATT",
            method = "nearest", distance = "mahalanobis")

のように入力します。method = "nearest"は最近傍マッチングを意味し、distance = "mahalanobis"はマハラノビス距離を意味します。estimand = "ATT"はATTを推定することを意味します。{MatchIt}の最近傍マッチングの場合、"ATT"、または"ATC"のみ指定可能です。早速やってみましょう。

MH_Matching <- matchit(treat ~ age + educ + black + hispanic + married + 
                           nodegree + re74 + re75, data = lalonde_df,
                       method = "nearest", distance = "mahalanobis",
                       estimand = "ATT")

バランスチェックは数値的に確認する方法と視覚的にする方法があります。まずは、数値的に確認してみましょう。

# バランスチェック (Numerical)
summary(MH_Matching)
## 
## Call:
## matchit(formula = treat ~ age + educ + black + hispanic + married + 
##     nodegree + re74 + re75, data = lalonde_df, method = "nearest", 
##     distance = "mahalanobis", estimand = "ATT")
## 
## Summary of Balance for All Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age            25.8162       28.0303         -0.3094     0.4400    0.0813
## educ           10.3459       10.2354          0.0550     0.4959    0.0347
## black           0.8432        0.2028          1.7615          .    0.6404
## hispanic        0.0595        0.1422         -0.3498          .    0.0827
## married         0.1892        0.5128         -0.8263          .    0.3236
## nodegree        0.7081        0.5967          0.2450          .    0.1114
## re74         2095.5737     5619.2365         -0.7211     0.5181    0.2248
## re75         1532.0553     2466.4844         -0.2903     0.9563    0.1342
##          eCDF Max
## age        0.1577
## educ       0.1114
## black      0.6404
## hispanic   0.0827
## married    0.3236
## nodegree   0.1114
## re74       0.4470
## re75       0.2876
## 
## 
## Summary of Balance for Matched Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age            25.8162       25.5351          0.0393     0.7451    0.0341
## educ           10.3459       10.4703         -0.0618     0.8918    0.0122
## black           0.8432        0.2595          1.6057          .    0.5838
## hispanic        0.0595        0.1297         -0.2971          .    0.0703
## married         0.1892        0.4486         -0.6625          .    0.2595
## nodegree        0.7081        0.6108          0.2140          .    0.0973
## re74         2095.5737     2545.7229         -0.0921     1.2602    0.0748
## re75         1532.0553     1645.4710         -0.0352     1.3141    0.0386
##          eCDF Max Std. Pair Dist.
## age        0.1243          0.3188
## educ       0.0973          0.2285
## black      0.5838          1.6949
## hispanic   0.0703          0.6172
## married    0.2595          1.0765
## nodegree   0.0973          0.2378
## re74       0.3405          0.2010
## re75       0.1838          0.1586
## 
## Percent Balance Improvement:
##          Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
## age                 87.3       64.2      58.1     21.2
## educ               -12.5       83.7      64.8     12.6
## black                8.8          .       8.8      8.8
## hispanic            15.1          .      15.1     15.1
## married             19.8          .      19.8     19.8
## nodegree            12.6          .      12.6     12.6
## re74                87.2       64.8      66.7     23.8
## re75                87.9     -511.2      71.2     36.1
## 
## Sample Sizes:
##           Control Treated
## All           429     185
## Matched       185     185
## Unmatched     244       0
## Discarded       0       0

表が3つでますが、上から(1)マッチング前の各共変量の情報、(2)マッチング後の各共変量の情報、(3)バランス改善の指標となります。(3)から見る場合、hispanic変数は改善率が100%となっています。実際、マッチング前は統制群と処置群におけるhispanic変数の平均値の差分は-0.0827でしたが、マッチング後は0となり、完全にマッチングされたことが分かります。他にも多くの共変量においてバランスの改善が見られます。eudcはむしろ悪化しましたが、他の変数はかなり改善されているので、このままいきましょう。

つづいて視覚的な方法です。

# バランスチェック (Graphical): QQプロットによるバランスチェック
plot(MH_Matching)

各ポイントが45度直線上に並ぶのが理想的なバランシングです。バランスが改善されていることが確認できます。

もっと直感的なバランスチェックの方法として{BalanceR}パッケージと同様、標準化差分を使う方法です。マッチング前後のバランスを同時に確認するには{cobalt}パッケージのlove.plotが便利です。

love.plot(MH_Matching, thresholds = 0.25, abs = TRUE)

thresholds引数は垂直線（破線）の位置、absは標準化差分を絶対値で示すことを意味します。マッチング後の標準化差分（Adjusted; 赤い点）が0.25より左側に位置している場合、バランスしていると判断できます2。他にもマッチング後の標準化差分がマッチング前（Unadjusted; 青い点）より改善されるいるか否かも判断できます。今回の例だと、大幅にバランスが改善されました。blackはまだバランスが取れておりませんが、それでも大幅に改善されていることが分かります。

それではATTを推定してみましょう。推定方法としてはノンパラメトリックな方法とパラメトリック方法がありますが、結果は変わりません。ノンパラメトリックな方法はペアごとの差分を計算し、その平均値を求める方法ですが、マッチング済みのデータに対し、処置変数を結果変数を回帰させることも、結果的には同じことを行うことになります。したがって、もっと簡単なパラメトリック方法でATTを推定します。ノンパラメトリック方法によるATTの推定は本資料の付録を参照してください。

回帰分析を行うためにはデータが必要ですね。つまり、マッチングされないケースをデータから除去する必要があります。したがって、マッチングされた後のデータセットを抽出する必要があります。使う関数はmatch.data()関数です。

# マッチング後データを抽出
MH_Data <- match.data(MH_Matching)

マッチングデータが取れたら、その中身を確認してみましょう。

dim(MH_Data)
## [1] 370  14
MH_Data

データのサイズは370行14列です。この370行には意味があります。それは処置群の大きさの2倍という点です。多くの場合、マッチングから計算される処置効果はATEではなく、ATTです。したがって、処置群のデータを100%活用し、統制群から共変量が最も近いケースを抽出&マッチングすることになります。だから、マッチング後のサンプルサイズは処置群のサイズの2倍になります。

それでは職業訓練のATTを推定します。方法は簡単です。このマッチング後のデータを用い、単回帰分析を行うのみです。

# 処置効果の確認 (Model-Based)
MH_Fit <- lm(re78 ~ treat, data = MH_Data)

modelsummary(MH_Fit)
Model 1
(Intercept) 5842.100
(528.119)
treat 507.043
(746.873)
Num.Obs. 370
R2 0.001
R2 Adj. -0.001
AIC 7624.8
BIC 7636.6
Log.Lik. -3809.419
F 0.461

処置効果は約507.043ドルであり、付録にあるノンパラメトリック推定値と一致していることが分かります。今回の結果は重回帰分析よりも推定値が低めであり、統計的に有意に職業訓練の効果があったとは言えないという結果が得られましたね。

ちなみに、{MatchIt}パッケージを使った最近傍マッチングのの結果は行う度に変化することがあります。なぜなら、{MatchIt}パッケージを使った最近傍マッチングの場合、処置群 (統制群)から一つのケースを選択し、最も近い統制群 (処置群)とマッチングさせます。マッチングされたケースは次のステップからはマッチング対象から除外されます3。また、1:1マッチングの場合4、同距離に複数のマッチング対象があると、ランダムに1つのみを選択します。最近傍マッチングを用いる際は、複数推定を行い、推定が安定するかを確認し、不安定な場合は他の手法を使うか、k-最近傍マッチングなどを使ってみましょう。

Coarsened Exact Matching

CEMはマハラノビス最近傍マッチング同様、matchit()関数を使います。ただし、cemパッケージをインストールしておく必要があります。まず、Exact Matching類の手法は距離を図る必要がないので、distance引数は不要です。method = ...引数を"nearest" (最近傍)から"cem"に替えるだけです。推定可能な処置効果はATE、ATT、ATCですが、ここではATTを推定してみましょう。

マッチングをしたらmatch.data()でマッチングされたデータを抽出します。

CEM_Matching <- matchit(treat ~ age + educ + black + hispanic + married + 
                            nodegree + re74 + re75, data = lalonde_df,
                       method = "cem", estimand = "ATT")

summary(CEM_Matching)
## 
## Call:
## matchit(formula = treat ~ age + educ + black + hispanic + married + 
##     nodegree + re74 + re75, data = lalonde_df, method = "cem", 
##     estimand = "ATT")
## 
## Summary of Balance for All Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age            25.8162       28.0303         -0.3094     0.4400    0.0813
## educ           10.3459       10.2354          0.0550     0.4959    0.0347
## black           0.8432        0.2028          1.7615          .    0.6404
## hispanic        0.0595        0.1422         -0.3498          .    0.0827
## married         0.1892        0.5128         -0.8263          .    0.3236
## nodegree        0.7081        0.5967          0.2450          .    0.1114
## re74         2095.5737     5619.2365         -0.7211     0.5181    0.2248
## re75         1532.0553     2466.4844         -0.2903     0.9563    0.1342
##          eCDF Max
## age        0.1577
## educ       0.1114
## black      0.6404
## hispanic   0.0827
## married    0.3236
## nodegree   0.1114
## re74       0.4470
## re75       0.2876
## 
## 
## Summary of Balance for Matched Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age            21.6462       21.2936          0.0493     0.8909    0.0130
## educ           10.3692       10.2795          0.0446     0.8894    0.0100
## black           0.8615        0.8615          0.0000          .    0.0000
## hispanic        0.0308        0.0308          0.0000          .    0.0000
## married         0.0308        0.0308          0.0000          .    0.0000
## nodegree        0.6769        0.6769          0.0000          .    0.0000
## re74          489.0612      697.6115         -0.0427     1.3916    0.0435
## re75          362.4365      520.8378         -0.0492     0.8501    0.0403
##          eCDF Max Std. Pair Dist.
## age        0.0897          0.1373
## educ       0.1000          0.1949
## black      0.0000          0.0000
## hispanic   0.0000          0.0000
## married    0.0000          0.0000
## nodegree   0.0000          0.0000
## re74       0.3118          0.0968
## re75       0.1598          0.1635
## 
## Percent Balance Improvement:
##          Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
## age                 84.1       85.9      84.0     43.1
## educ                18.8       83.3      71.2     10.2
## black              100.0          .     100.0    100.0
## hispanic           100.0          .     100.0    100.0
## married            100.0          .     100.0    100.0
## nodegree           100.0          .     100.0    100.0
## re74                94.1       49.7      80.6     30.2
## re75                83.0     -263.5      70.0     44.4
## 
## Sample Sizes:
##               Control Treated
## All            429.       185
## Matched (ESS)   41.29      65
## Matched         75.        65
## Unmatched      354.       120
## Discarded        0.         0

CEMの場合、マッチングされないブロックは捨てられるため、マハラノビス距離最近傍マッチングよりもサンプルサイズが小さくなります。また、処置群と統制群のサンプルサイズも不均衡になります。マッチング結果を見ると、処置群からは65ケース、統制群からは75サンプルのみ残ることになりました。

つづいて、cobalt::love.plot()を使用して、バランスチェックを行います。

love.plot(CEM_Matching, thresholds = 0.25, abs = TRUE)

CEMの場合、マハラノビス距離最近傍マッチングよりもバランスが大きく改善されたことが分かります。その理由は簡単です。最近傍マッチングの場合、最も近いケースであれば、どれほど離れていてもマッチングされます。一方、CEMは正確マッチングの一種であるため、ある程度離れているケースを捨ててしまうからです。結局は共変量が非常に近いケースのみを残すことになります。

それでは、match.data()関数を使ってマッチング後のデータを抽出してみましょう。

CEM_Data <- match.data(CEM_Matching)

処置効果の推定の前に、抽出されたデータの中身を見てましょう。元のデータにはweightsという列はありませんでしたが、いきなりできました。実はこれが、CEMに使う重み変数です。自動的に計算してデータセットに追加してくれたわけです5。

CEM_Data

処置効果の推定は重み付き回帰分析 (WLS)で行います。既存のlm()関数にweight = ...引数を追加し、データ内の重み列をしてするだけです。

# 処置効果の確認
# weight = Wが追加されていることを確認
CEM_Fit <- lm(re78 ~ treat, data = CEM_Data, weights = weights)

# 処置効果の確認
modelsummary(CEM_Fit)
Model 1
(Intercept) 5265.785
(850.457)
treat 1070.907
(1248.129)
Num.Obs. 140
R2 0.005
R2 Adj. -0.002
AIC 2924.7
BIC 2933.5
Log.Lik. -1459.329
F 0.736

今回、処置効果の推定値は約1070.907ドルです。

傾向スコアマッチング

傾向スコアマッチングも、これまでのコマンドとほぼ同じです。マハラノビス最近傍マッチングのコマンドからdistance = ...引数を抜けば、傾向スコアマッチングができます6。また、replace引数にTRUEを指定します。これは一度マッチングしたケースを捨てずに、次のマッチングにも使うことを意味します。デフォルトはFALSEとなっており、この場合、一度マッチングされたケースは二度とマッチングされません。したがって、replace = TRUEを指定した方が共変量バランスがより改善されます。ただし、有効サンプルサイズ（Effective Sample Size; ESS）が小さくなり、精度が悪くなる点、場合によっては特殊な標準誤差7を使う必要があるといった欠点があります。これは最近傍マッチングに共通する問題ですが、ここではこれまでの結果と比較するために既定値のままにしましょう。

PS_Matching <- matchit(treat ~ age + educ + black + hispanic + married + 
                           nodegree + re74 + re75, data = lalonde_df,
                       method = "nearest", estimand = "ATT")

summary(PS_Matching)
## 
## Call:
## matchit(formula = treat ~ age + educ + black + hispanic + married + 
##     nodegree + re74 + re75, data = lalonde_df, method = "nearest", 
##     estimand = "ATT")
## 
## Summary of Balance for All Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance        0.5774        0.1822          1.7941     0.9211    0.3774
## age            25.8162       28.0303         -0.3094     0.4400    0.0813
## educ           10.3459       10.2354          0.0550     0.4959    0.0347
## black           0.8432        0.2028          1.7615          .    0.6404
## hispanic        0.0595        0.1422         -0.3498          .    0.0827
## married         0.1892        0.5128         -0.8263          .    0.3236
## nodegree        0.7081        0.5967          0.2450          .    0.1114
## re74         2095.5737     5619.2365         -0.7211     0.5181    0.2248
## re75         1532.0553     2466.4844         -0.2903     0.9563    0.1342
##          eCDF Max
## distance   0.6444
## age        0.1577
## educ       0.1114
## black      0.6404
## hispanic   0.0827
## married    0.3236
## nodegree   0.1114
## re74       0.4470
## re75       0.2876
## 
## 
## Summary of Balance for Matched Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance        0.5774        0.3629          0.9739     0.7566    0.1321
## age            25.8162       25.3027          0.0718     0.4568    0.0847
## educ           10.3459       10.6054         -0.1290     0.5721    0.0239
## black           0.8432        0.4703          1.0259          .    0.3730
## hispanic        0.0595        0.2162         -0.6629          .    0.1568
## married         0.1892        0.2108         -0.0552          .    0.0216
## nodegree        0.7081        0.6378          0.1546          .    0.0703
## re74         2095.5737     2342.1076         -0.0505     1.3289    0.0469
## re75         1532.0553     1614.7451         -0.0257     1.4956    0.0452
##          eCDF Max Std. Pair Dist.
## distance   0.4216          0.9740
## age        0.2541          1.3938
## educ       0.0757          1.2474
## black      0.3730          1.0259
## hispanic   0.1568          1.0743
## married    0.0216          0.8281
## nodegree   0.0703          1.0106
## re74       0.2757          0.7965
## re75       0.2054          0.7381
## 
## Percent Balance Improvement:
##          Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
## distance            45.7     -239.6      65.0     34.6
## age                 76.8        4.6      -4.2    -61.1
## educ              -134.8       20.4      31.2     32.1
## black               41.8          .      41.8     41.8
## hispanic           -89.5          .     -89.5    -89.5
## married             93.3          .      93.3     93.3
## nodegree            36.9          .      36.9     36.9
## re74                93.0       56.8      79.1     38.3
## re75                91.2     -800.7      66.3     28.6
## 
## Sample Sizes:
##           Control Treated
## All           429     185
## Matched       185     185
## Unmatched     244       0
## Discarded       0       0

処置効果の推定の前に、バランスが取れているかどうかを確認してみましょう。これまではQQプロットのみでバランスチェックをしましたが、今回はヒストグラムとJitterプロットで確認します。

# 1. QQプロットによるバランスチェック
plot(PS_Matching)

QQプロットを見ると、バランスが多少改善されていることが分かります。

# 2. ヒストグラムによるバランスチェック
plot(PS_Matching, type = "hist")

ヒストグラムからもバランスが改善されたことが確認できますね。出来れば、ヒストグラムを重ねて比較したいんですが、残念です。あとでオーバーラップされたヒストグラムの作成方法について解説します。

# 3. 散布図によるバランスチェック
par(mfrow = c(1, 1))
plot(PS_Matching, type = "jitter")

## [1] "To identify the units, use first mouse button; to stop, use second."
## integer(0)

散布図をみると、傾向スコア (\(e\))が小さいケースがマッチングから除外されたことが分かります。実際、処置群には傾向スコアが小さいケースが非常に少ないため、これは妥当な結果でしょう。

むろん、おなじみのcobalt::love.plot()でもバランスチェックができます。

love.plot(PS_Matching, thresholds = 0.25, abs = TRUE)

それでは、ATT推定のためにマッチング後のデータを抽出します。

# 傾向スコアマッチング後のデータセットを抽出
PS_Data <- match.data(PS_Matching)

傾向スコアを用いたATTをの推定もこれまでと同様、回帰分析を用います。普通の回帰分析ですが、用いるデータは傾向スコアでマッチングされたデータのみになります8。

# 処置効果の推定
PS_Fit <- lm(re78 ~ treat, data = PS_Data)

modelsummary(PS_Fit)
Model 1
(Intercept) 5454.776
(516.396)
treat 894.367
(730.295)
Num.Obs. 370
R2 0.004
R2 Adj. 0.001
AIC 7608.2
BIC 7620.0
Log.Lik. -3801.114
F 1.500

処置効果は約894.367ドルでした。

ちなみに、講義スライド47〜48ページ目の例を{MatchIt}パッケージを使うと以下のようなコードになります。講義スライドでは一度マッチングされたケースを再利用したため、replace = TRUEを指定し、ATT/ATCを推定する際にweights =引数を指定する必要があります。

# 講義スライドのデータ
PS_df2 <- read.csv("Data/Day2_Data2.csv")

p47_mat <- matchit(Treat ~ Age + Edu, estimand = "ATT",
                   data = PS_df2, method = "nearest", replace = TRUE)
p48_mat <- matchit(Treat ~ Age + Edu, estimand = "ATC",
                   data = PS_df2, method = "nearest", replace = TRUE)

p47_data <- match.data(p47_mat)
p48_data <- match.data(p48_mat)

p47_fit <- lm(Outcome ~ Treat, data = p47_data, weights = weights)
p48_fit <- lm(Outcome ~ Treat, data = p48_data, weights = weights)

modelsummary(list("ATT" = p47_fit, "ATC" = p48_fit))
ATT ATC
(Intercept) 4.909 4.308
(0.830) (0.599)
Treat 1.818 3.923
(1.031) (1.012)
Num.Obs. 17 20
R2 0.172 0.455
R2 Adj. 0.116 0.425
AIC 76.7 93.8
BIC 79.2 96.8
Log.Lik. -35.343 -43.906
F 3.108 15.031

手計算から得られた結果と一致することが分かります。

IPW推定量

最後にATTでなく、ATEが推定の対象となるIPW推定量を計算してみましょう。{WeightIt}パッケージを使用すれば便利ですが、ここではあえて使わずにIPW推定量を計算してみましょう。こうすることでIPWが具体的にどのような仕組みで推定されるかが理解できるはずです。一通りの推定が終わりましたら{WeightIt}パッケージの使い方も紹介します。

パッケージを使わない方法

まずは、傾向スコアの算出ですが、ここではロジスティック回帰分析で傾向スコアを計算します9。

# Logitで傾向スコアを推定
PS <- glm(treat ~ age + educ + black + hispanic + married + nodegree +
              re74 + re75, data = lalonde_df,
          family = binomial("logit"))

続いて、既存のデータフレームに傾向スコアが格納された列を追加します。まず、既存のlalondeデータセットをMatching.dfという名前で複製します。

新しいデータフレームにPSという列を作成し、ここに傾向スコアを格納します。傾向スコアはロジスティック回帰分析のオブジェクトに$fitted.valuesを加えるだけで抽出できます。

# lalondeデータフレームを複製
ipw_df <- lalonde_df

# 傾向スコア推定値をデータフレーム内に格納
ipw_df$PS <- PS$fitted.values

次は、各ケースの重みを計算し、Wという新しい列として追加します。傾向スコアは、もしケースが処置群 (treat == 1)なら、\(\frac{1}{e}\)を、統制群 (treat == 0)なら\(\frac{1}{(1-e)}\)となります。

# 重み変数を作成
## ケースが処置群なら、重み = 傾向スコアの逆数
## ケースが統制群あら、重み = (1 - 傾向スコア)の逆数
ipw_df <- ipw_df %>%
    mutate(W = ifelse(treat == 1, (1 / PS), (1 / (1 - PS))))

以上のコードは以下のコードでも同じく作動します。

ipw_df <- ipw_df %>%
    mutate(W = treat * (1 / PS) + (1 - treat) * (1 / (1 - PS)))

むろん、dplyrのmutate()を使わずに、以下のように入力しても重み変数は作成可能である。

ipw_df$W <- ifelse(ipw_df$treat == 1,
                   1 / ipw_df$PS,
                   1 / (1 - ipw_df$PS))

最後に処置効果の推定です。CEMの時と同じやり方です。weight = ...引数を追加するだけです。

# 回帰分析にweight = ...を指定するだけ
IPW_fit <- lm(re78 ~ treat, data = ipw_df, weights = W)

modelsummary(IPW_fit)
Model 1
(Intercept) 6422.839
(397.427)
treat 224.676
(577.657)
Num.Obs. 614
R2 0.000
R2 Adj. -0.001
AIC 12796.0
BIC 12809.3
Log.Lik. -6394.998
F 0.151

推定された処置効果は約224.676ドルです。

{WeightIt}パッケージを使用した場合

それでは{WeightIt}パッケージを使ってみましょう。これまで使ってきました{MatchIt}パッケージや傾向スコア算出のためのglm()関数と使い方が非常に似ています。まず、第一引数として処置変数を結果変数、処置有無に影響を与えると考えられる共変量を説明変数とした数式オブジェクト（formulaクラス）を指定します。続いて、データ（data）、IPW算出の方法（method）、推定の対象（estimand）を指定します。データはdata = lalonde_dfとし、傾向スコアからIPWを算出するためmethod = "ps"を指定、最後にATT推定のためにestimand = "ATT"を指定します。ATEあるいはATCを計算する場合は"ATT"を適宜変更してください。また、IPW算出のために今回は傾向スコアを使いましたが、「処置を受ける確率」が計算できるなら何でも良いです。たとえば、Imai and Ratkovic (2014)10が推奨しているCovariate Balancing Propensity Score (CBPS) を使用する場合は"ps"の代わりに"cbps"を、複数の推定を組み合わせるスーパーラーニングをする場合は"super"11などが使えます。他にもエントロピーバランシングなど様々なオプションが提供されています。ここでは、先程の手法との比較のために、デフォルトの"ps"を使います。

Weighted_Data <- weightit(treat ~ age + educ + black + hispanic + married +
                            nodegree + re74 + re75, data = lalonde_df,
                          method = "ps", estimand = "ATT")

{WeightIt}パッケージの便利なところはcobalt::love.plot()によるバランスチェックが簡単ということです。実際にやってみましょう。

love.plot(Weighted_Data, thresholds = 0.25, abs = TRUE)

最後にIPW推定量を計算してみましょう。ここは先ほどの面倒くさい方法と同じですが、重み付け変数はweightit()関数から得られたオブジェクト（Weighted_Data）のweights列を使用します。推定ができたら結果を出力します。

IPW_Result <- lm(re78 ~ treat, data = lalonde_df,
                 weights = Weighted_Data$weights)

modelsummary(IPW_Result)
Model 1
(Intercept) 5135.072
(401.194)
treat 1214.071
(568.904)
Num.Obs. 614
R2 0.007
R2 Adj. 0.006
AIC 13241.9
BIC 13255.1
Log.Lik. -6617.942
F 4.554

という結果が得られましたが、これは{WeightIt}パッケージを使わなかった結果と完全に一致しますね。

付録

ノンパラメトリック推定

ここではマハラノビス最近傍マッチングの例を使って、ATTを推定してみましょう。まずは、処置群の潜在的結果を統制群から割り当てて差分を算出し、その平均値を求るというのがノンパラメトリックな推定方法です。以下のコードが理解できるようになれば、他の最近傍マッチング（傾向スコアマッチングなど）も外部パッケージを使わずに計算できるようになります。

各処置群がどのケースとマッチングされたかはマッチングのオブジェクト名の後に$match.matrixを付けることで抽出できます。ちょっとやってみましょう。

# 全部表示させると長くなんるので、最初の6行まで確認
head(MH_Matching$match.matrix) 
##   [,1] 
## 1 "553"
## 2 "557"
## 3 "478"
## 4 "567"
## 5 "579"
## 6 "438"

左が処置群、右が統制群ですね。それぞれケースのIDを抽出しましょう。

# MH.Matching$match.matrixは1列の行列です。
# 処置群IDは行の名前として格納され、統制群IDは1列目に入っています。
MH_Treat   <- as.vector(row.names(MH_Matching$match.matrix))
MH_Control <- as.vector(MH_Matching$match.matrix[, 1])

# 中身を確認
print(MH_Treat)
##   [1] "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "9"   "10"  "11"  "12" 
##  [13] "13"  "14"  "15"  "16"  "17"  "18"  "19"  "20"  "21"  "22"  "23"  "24" 
##  [25] "25"  "26"  "27"  "28"  "29"  "30"  "31"  "32"  "33"  "34"  "35"  "36" 
##  [37] "37"  "38"  "39"  "40"  "41"  "42"  "43"  "44"  "45"  "46"  "47"  "48" 
##  [49] "49"  "50"  "51"  "52"  "53"  "54"  "55"  "56"  "57"  "58"  "59"  "60" 
##  [61] "61"  "62"  "63"  "64"  "65"  "66"  "67"  "68"  "69"  "70"  "71"  "72" 
##  [73] "73"  "74"  "75"  "76"  "77"  "78"  "79"  "80"  "81"  "82"  "83"  "84" 
##  [85] "85"  "86"  "87"  "88"  "89"  "90"  "91"  "92"  "93"  "94"  "95"  "96" 
##  [97] "97"  "98"  "99"  "100" "101" "102" "103" "104" "105" "106" "107" "108"
## [109] "109" "110" "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
## [121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130" "131" "132"
## [133] "133" "134" "135" "136" "137" "138" "139" "140" "141" "142" "143" "144"
## [145] "145" "146" "147" "148" "149" "150" "151" "152" "153" "154" "155" "156"
## [157] "157" "158" "159" "160" "161" "162" "163" "164" "165" "166" "167" "168"
## [169] "169" "170" "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
## [181] "181" "182" "183" "184" "185"
print(MH_Control)
##   [1] "553" "557" "478" "567" "579" "438" "592" "528" "590" "564" "540" "566"
##  [13] "561" "608" "463" "552" "556" "544" "578" "547" "460" "473" "585" "537"
##  [25] "596" "563" "603" "524" "577" "525" "586" "604" "516" "558" "454" "581"
##  [37] "523" "571" "457" "527" "568" "422" "583" "450" "573" "359" "449" "451"
##  [49] "569" "570" "532" "565" "559" "374" "601" "594" "411" "600" "575" "538"
##  [61] "550" "420" "382" "384" "613" "468" "421" "464" "539" "584" "477" "588"
##  [73] "599" "430" "415" "396" "533" "526" "409" "574" "377" "426" "546" "609"
##  [85] "447" "529" "397" "508" "406" "509" "506" "598" "610" "522" "560" "513"
##  [97] "555" "593" "365" "429" "462" "458" "459" "542" "424" "511" "452" "562"
## [109] "400" "474" "504" "576" "520" "456" "445" "507" "344" "391" "443" "366"
## [121] "405" "461" "387" "497" "475" "351" "501" "370" "363" "399" "436" "311"
## [133] "362" "317" "352" "340" "410" "390" "353" "280" "442" "302" "333" "304"
## [145] "329" "321" "325" "338" "303" "319" "275" "301" "327" "328" "323" "339"
## [157] "343" "314" "259" "265" "320" "324" "277" "291" "257" "232" "271" "255"
## [169] "251" "247" "274" "261" "217" "235" "243" "220" "233" "269" "216" "202"
## [181] "192" "208" "191" "186" "193"

次は、元のlalondeデータから処置群と統制群の結果変数 (re78)を抽出し、MH.Matched.dfというデータフレームとしてまとめます。続いて、各ケースのITEを計算します。

MH_Matched_df <- tibble(
    Treat_ID   = MH_Treat,
    Control_ID = MH_Control,
    Treat_Y    = lalonde[MH_Treat,]$re78,
    Control_Y  = lalonde[MH_Control,]$re78
)

MH_Matched_df <- MH_Matched_df %>%
    mutate(ITE = Treat_Y - Control_Y)
MH_Matched_df

1行目を見ると1は553とマッチングされ、それぞれの結果変数は9930.046ドル、0ドルです。その差分がITEであり、9930.046ドルですね。

これらITEの平均値がATTとなります。早速計算してみましょう。

MH_ATT <- mean(MH_Matched_df$ITE)

MH_ATT
## [1] 507.0434

ATTは約507.043ドルでした。回帰分析でやった結果と一致しますね。大人しく回帰分析でやりましょう。

バランスチェック

これまでは主にQ-Q Plotを用いたバランスチェックを紹介しましたが、MatchItパッケージは標準化差分/標準化バイアスのバランスチェックにも対応します。やり方としては、

X <- matchit(...)
X_Blc <- summary(X, standardize = TRUE)

plot(X, interactive = FALSE)

です。たとえば、傾向スコアマッチングの例だと、

PS_Blc <- summary(PS_Matching, standardize = TRUE)

plot(PS_Blc, interactive = FALSE)
abline(v = 0.25, col = "red")

interactive = TRUEにすると各ポイントがどの変数に対応するかが分かります。灰色のポイントはバランスが改善された共変量であり、黒は悪化した共変量です。また、これは前の講義で解説した標準化差分ですが、100を掛けていないものであり、なお、絶対値です。したがって、バランスの基準は0.1となります。多くの共変量においてバランスが改善が見られますが、0.1未満の共変量も少なく、大きな改善とは考えられませんね。

もっと分かりやすい図がほしければ、これまで何回か使ってきましたcobaltパッケージがおすすめです。

love.plot(PS_Matching, threshold = 0.25, abs = TRUE)

赤い点がマッチング前の標準化差分、青い点がマッチング後の標準化差分です。図の上にあるdistanceは傾向スコアを意味します。educとhispanはマッチングによってバランスが悪化しましたが、他の変数が改善された幅を考えると許容範囲かもしれません。それでも、0.25の垂直線の右側にある青い点が3つもあり、傾向スコアの標準化差分もかなり大きいのが分かります。

他にも傾向スコアをマッチング前後に分けてヒストグラム化することもいいでしょう。

bal.plot(PS_Matching, var.name = "distance", which = "both",
         type = "histogram", mirror = TRUE)

左がマッチング前、右が後です。また、赤いヒストグラムは統制群、青は処置群です。もし、傾向スコアのバランスが取れているならヒストグラムは上下に対称となります。マッチング後の傾向スコアの標準化差分はかなり大きかったのですが、それでも大幅の改善は見られました。少なくとも、普通に回帰分析するよりも良さそうですね。

複数の手法を用いた場合、推定値の比較

これまでの推定値を同時に示す図を作ってみましょう。そのためには各推定結果のオブジェクトから処置変数の係数を抽出する必要があります。たとえば、重回帰分析によるATEの推定結果はRaw_Fit2に格納されています。lm()やglm()関数で行った推定結果から係数のみを抽出する関数はcoef()です。

coef(Raw_Fit2)
##   (Intercept)         treat           age          educ         black 
##    66.5145197  1548.2438020    12.9776337   403.9412309 -1240.6440820 
##      hispanic       married      nodegree          re74          re75 
##   498.8968650   406.6208454   259.8173680     0.2963774     0.2315259

各係数の推定値がベクトルとして表示されます。ここで処置変数（treat）の係数はベクトルの2番目の要素であるため、coef(Raw_Fit2)[2]で抽出できます。他の手法による推定についても同じです。なぜなら、これまで処置変数を1番目の説明変数として使ってきたからです。

一方、標準誤差の抽出はちょっと面倒くさいです。いくつかの方法がありますが、ここではsummary(オブジェクト名)$coefを使います。推定結果を出力するsummary()関数のcoef要素を抽出すると以下のようなデータフレームが抽出されます。

summary(Raw_Fit2)$coef
Estimate Std. Error t value Pr(>|t|)
(Intercept) 66.5145197 2436.7458028 0.0272965 0.9782323
treat 1548.2438020 781.2792975 1.9816778 0.0479682
age 12.9776337 32.4889061 0.3994482 0.6897042
educ 403.9412309 158.9062381 2.5420099 0.0112704
black -1240.6440820 768.7644071 -1.6138157 0.1070897
hispanic 498.8968650 941.9425489 0.5296468 0.5965515
married 406.6208454 695.4723228 0.5846686 0.5589889
nodegree 259.8173680 847.4420524 0.3065901 0.7592610
re74 0.2963774 0.0582726 5.0860483 0.0000005
re75 0.2315259 0.1046199 2.2130198 0.0272694

ここで、treatの標準誤差は2行2列目に格納されていることが分かります。したがって、treat標準誤差はsummary(Raw_Fit2)$coef[2, 2]で抽出できます。他にもsummary(Raw_Fit2)から抽出できる要素は色々あります。興味のある方はstr(summary(Raw_Fit2))で確認してみてください。

Comparison_df <- data.frame(
    Method = c("Single Regression",
               "Multiple Regression",
               "Mahalanobis (ATT)",
               "CEM (ATT)",
               "Propensity Score (ATT)",
               "IPW (ATT)"),
    Est    = c(coef(Raw_Fit1)[2],
               coef(Raw_Fit2)[2],
               coef(MH_Fit)[2],
               coef(CEM_Fit)[2],
               coef(PS_Fit)[2],
               coef(IPW_fit)[2]),
    SE     = c(summary(Raw_Fit1)$coef[2, 2],
               summary(Raw_Fit2)$coef[2, 2],
               summary(MH_Fit)$coef[2, 2],
               summary(CEM_Fit)$coef[2, 2],
               summary(PS_Fit)$coef[2, 2],
               summary(IPW_fit)$coef[2, 2])
)

続いて、各手法名（Method列）を順序付きの変数に変換します。これをしないと、各手法がアルファベット順で並ぶことになります。順番が、データフレームに登場する順番の逆となります。したがって、fct_inorder()を使い、「登場順でfactor化」をし、fct_rev()で順番の反転を行います。fct_*()関数群はforcatsパッケージが提供しており、詳細は宋・矢内の第14章を参照してください。

Comparison_df <- Comparison_df %>%
  mutate(Method = fct_inorder(Method),
         Method = fct_rev(Method))

続きまして、95%信頼区間の下限（LL）と上限（UL）を追加し、\(\alpha = 0.05\)水準で統計的有意か否かを示す変数も作成します。

Comparison_df <- Comparison_df %>%
    mutate(LL  = Est + qnorm(0.025) * SE,
           UL  = Est + qnorm(0.975) * SE,
           Sig = if_else(LL * UL > 0 , "Significant", "Insignificant"))
Comparison_df
Method Est SE LL UL Sig
Single Regression -635.0262 657.1374 -1922.99190 652.9395 Insignificant
Multiple Regression 1548.2438 781.2793 16.96452 3079.5231 Significant
Mahalanobis (ATT) 507.0434 746.8731 -956.80109 1970.8878 Insignificant
CEM (ATT) 1070.9073 1248.1293 -1375.38120 3517.1959 Insignificant
Propensity Score (ATT) 894.3675 730.2948 -536.98398 2325.7189 Insignificant
IPW (ATT) 224.6763 577.6575 -907.51150 1356.8641 Insignificant

宋のように「表より図がいい！」という方は以下のように可視化することも可能です。

Comparison_df %>%
    ggplot() +
    # x = 0の垂直線を引く。線の種類は破線 (linetype = 2)で
    geom_vline(xintercept = 0, linetype = 2) +
    # 点と線を同時に出力するgeom_pointrange幾何オブジェクトを使用
    # 点の位置 (x, y)はEst, Methodに、線の下限と上限 (xmin, xmax)はLL, UL
    # Sigの値によって色分けするためにcolor引数もマッピング
    geom_pointrange(aes(x = Est, y = Method, xmin = LL, xmax = UL,
                        color = Sig), size = 0.75) +
    # colorにマッピングされているSigの値がSignificantなら色をblackに
    # Insignificantならgray70にする
    scale_color_manual(values = c("Significant"   = "black",
                                  "Insignificant" = "gray70")) +
    # ラベルを付ける
    labs(y = "Methods", x = "Estimates (ATE or ATT)",
         color = "") +
    # minimalテーマを適用
    theme_minimal() +
    # 文字の大きさを調整
    theme(text = element_text(size = 16))

  1. lalondeデータセットを提供するパッケージは複数あり、それぞれデータセットの構成に違いがあります。どれを使っても実習には問題ございませんが、以下の内容を再現される場合はdata()内にpackage = "cobalt"を指定するようにしましょう。↩︎

  2. ここでは基準としている標準化差分は0.25ですが、{BalanceR}の25と同じです。↩︎

  3. これはmatchit()内にreplacement = TRUEを指定することで防ぐことができます。既定値はFALSEですが、TRUEを指定すれば、マッチングされた統制群ケースを2回以上マッチングすることができます。↩︎

  4. これはmatchit()内にratio引数を指定することで防ぐことができます。既定値は1であり、この場合は1:1マッチングです。↩︎

  5. これは重み変数を使わない最近傍マッチングをしても追加されます。ただし、その場合、weightsは全て1になります。↩︎

  6. 実はnearest = "logit"が省略されています。つまり、ロジスティック回帰分析から得られた傾向スコアの距離に基づくマッチングを意味します。↩︎

  7. Austin, Peter C., and Guy Cafri. 2020. “Variance Estimation When Using Propensity-Score Matching with Replacement with Survival or Time-to-Event Outcomes.” Statistics in Medicine, 39 (11): 1623–40.↩︎

  8. ただし、replacement = TRUEを指定した場合、重み付け回帰分析をする必要があります。。なぜなら、統制群から1つのケースが何回もマッチングされる可能性があるからです。統制群のあるケースが数回マッチングされてもPS_Dataにそのケースは1つしかありません。この場合は重みをより高く付ける必要があります。幸い、{MatchIt}パッケージはその重みの自動的に計算してくれます。具体的にはCEMと同様、weights変数がPS_Dataに追加されます。↩︎

  9. ロジスティック回帰分析ではなくても問題ありません。他、プロッビット、サポートベクトルマシン、回帰木、判別分析、ニューラル・ネットワークなどもあり、これらを組わせることも可能です。↩︎

  10. Imai, Kosuke and Marc Ratkovic. 2014. “Covariate Balancing Propensity Score.” Journal of the Royal Statistical Society, Series B, Vol. 76, No. 1, pp. 243-246.↩︎

  11. Super Learnerを使った例は「Cyrus, Samii, Laura Paler, and Sarah Zukerman Daly. 2016. “Retrospective Causal Inference with Machine Learning Ensembles: An Application to Anti-recidivism Policies in Colombia.” Political Analysis, 22 (4) pp. 434-456」を、日本語による解説は誰かの報告スライドを参照して下さい。↩︎